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We study the carrier-mediated exchange interaction, the so-called RKKY coupling, between two 
magnetic impurity moments in graphene using exact diagonalization on the honeycomb lattice. By 
using the tight-binding nearest neighbor band structure of graphene we also avoid the use of a 
momentum cut-off which plagues perturbative results in the Dirac continuum model formulation. 
We extract both the short and long impurity-impurity distance behavior and show on a qualitative 
agreement with earlier perturbative results in the long distance limit but also report on a few new 
findings. In the bulk the RKKY coupling is proportional to 1/|R|'^ and displays (1 -I- cos(2k_D • R.)- 
type oscillations. A-A sublattice coupling is always ferromagnetic whereas A-B subattice coupling 
is always antiferromagnetic and three times as large. We also study the effect of edges in zigzag 
graphene nanoribbons (ZGNRs). We find that for impurities on the edge the RKKY coupling 
decays exponentially because of the localized zero energy edge states and we also conclude that a 
non-perturbative treatment is essential for these edge impurities. For impurities inside a ZGNR the 
bulk characteristics are quickly regained. 

PACS numbers: 75.20.Hr, 75.75.-c, 73.20.-r 



I. INTRODUCTION 

Graphene is a two-dimensional honeycomb lattice of 
carbon and has since its isolation in 2004i. generated a lot 
of attention (see e.g. Ref. |3 and references therein). Its 
two-dimensionality, linear energy dispersion, where the 
quasiparticles are massless Dirac fermions, and chemi- 
cal potential tunable by a gate voltage are all novel fea- 
tures in an, essentially table-top, condensed matter sys- 
tem. Together with a very high mobility these prop- 
erties have helped to raise the expectation of graphene 
being a post-silicon era candidate.—"— For this, function- 
alization of graphene using for example finite geome- 
tries, adatoms, hydrogen chemisorption, or intrinsic de- 
fects such as vacancies, has become an important goal. 
Adding magnetic atoms or defects have the added benefit 
of opening the possibility for spintronics where not only 
the electron charge but also its spin is actively used in 
deviceS)^!^ One of the most important property of mag- 
netic impurities is the effective interaction between them 
propagated by the conduction electrons in the bulk host, 
the so-called Ruderman-Kittel-Kasuya-Yoshida (RKKY) 
couplingi^"— This coupling is crucial for magnetic order- 
ing of the impurities but also offers access to the intrinsic 
magnetic properties of the host. 

Previous work on the RKKY coupling in graphenei^r— 
have exclusively used a field-theory continuum model 
where the graphene band structure is approximated with 
a Dirac spectrum at each of the two inequivalent Bril- 
loiun zone corners. Using perturbation theory the RKKY 
coupling has then been calculated analytically from the 
static spin susceptibility of this model. The earliest 
worki2iH failed to recognize the importance of the bi- 
partite lattice and predicted that the RKKY coupling is 
always ferromagnetic. Later it was, however, shown by 
Saremil^ that for all bipartite lattices at half-filling the 
RKKY coupling is ferromagnetic (FM) only for impuri- 



ties on the same sublattice but antiferromagnetic (AFM) 
for impurities on different sublattices. A perturbative ap- 
proach also requires an explicit use of an ultraviolet mo- 
mentum cut-off scheme for the non-interacting graphene 
band structure and it was recently showed that a regular 
sharp cut-off does not produce the correct results thus 
demonstrating the need for a carefully chosen regulariza- 
tion schemCf^ This most likely explains the discrepancies 
in the later results for the RKKY couplingJ^ii^ Further- 
more, any results from a continuum model does not fully 
resolve the lattice structure which are likely to be im- 
portant for short impurity-impurity distances R. In fact, 
the later work-i^^iS, on the RKKY coupling in graphene 
have only considered the long-distance limit. While this 
is the traditional RKKY limit, knowledge of the short 
distance behavior is important in nano-structures as well 
as when comparing with ab-initio results where the unit 
cell always has a finite size. 

To circumvent the use of perturbation theory, the ul- 
traviolet cut-off dependency, and to also get results for 
finite impurity distances we will here explicitly calculate 
the RKKY coupling on the honeycomb lattice using exact 
diagonalization in a finite system. We show that by using 
a large enough system it is possible to extract both short- 
range and long-range behavior of the RKKY coupling. 
We report on a qualitative agreement between our results 
and one of the perturbative results in the long distance 
limifl-, but also report on a few new findings, includ- 
ing a phase shift for different sublattice coupling. This 
establishes not only that the standard perturbative ap- 
proach to RKKY coupling in graphene is in general valid 
but also finally settles the issue of the exact form of the 
RKKY coupling. Furthermore, we also study the RKKY 
coupling in zigzag graphene nanoribbons (ZGNRs) and 
show that the zigzag edge will significantly modify the re- 
sults due to the presence of the localized zero-energy edge 
statCFiir— In contrast to the bulk, a non-perturbative 
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treatment seems to be essential for impurities along the 
zigzag edge. For impurities inside a ZGNR bulk-like be- 
havior is however achieved even for narrow ribbons. 

The rest of the paper is organized as follows: In 
the next section we introduce the model and solution 
method. Then we discuss the results for on-site and pla- 
quette impurities, followed by the results for impurities 
on the edge and inside ZGNRs. We end with comparing 
our results with the previous perturbative results, and 
a discussion on how to experimentally achieve magnetic 
impurities in graphene as well as on the importance of 
electron-electron interactions . 



II. METHOD 

We are here focusing on how impurity magnetic mo- 
ments, or spins, interact with each other on a graphene 
surface. We are not concerned with the details of the in- 
teraction between the moments and graphene nor about 
how the moments are originally formed. We will there- 
fore model the interaction between an impurity moment 
and graphene with a simple Kondo coupling term. Using 
the nearest neighbor tight-binding Hamiltonian for the 
Pz-orbitals in graphene, the system can be formulated as 

H=-tJ2 Kabja + H.C.) + Jfc ^ S, • S„ (1) 

<i,j>,CT i— imp 

where aia- {bia) annihilates an electron on sublattice A 
(B) in unit cell i [see Fig. [IJa)], < i,j > means nearest 
neighbors, and a is the spin index. Moreover, S — ±Sz is 
the impurity spin and s = ^al^aapaii, with ctq,^ being the 
Pauli matrices, is the electron spin (for sublattice B in- 
terchange a ^ b). The constants entering are the nearest 
neighbor hopping in graphene t ~ 2.7 eV and the Kondo 
coupling Jk which depends on the particular impurity 
moment. We consider here only undoped graphene and 
thus no chemical potential term enters in Eq. p]). 




FIG. 1: (Color online) (a) The graphene honeycomb lattice 
with the two sublattices A and B in dark and light colors, re- 
spectively. The lattice unit vectors ci and C2 , lattice constant 
a — 2.46 A as well as the two most common directions, the 
zigzag and the armchair, are displayed, (b) Unit cell setup in 
the AFM configuration with the distance R between impurity 
spins shown as well as the padding p surrounding them. 

In standard RKKY perturbation theory^ the leading 



interaction between two impurity moments at sites i and 
j is given by 

Hrkky — JijSi ■ Sj, (2) 

with the effective RKKY coupling constant Jij propor- 
tional to the static spin susceptibility of the imbedding 
bulk, J^j oc Xij- 

A. Exact diagonalization 

Instead of using the above perturbative result for the 
RKKY coupling we will calculate Jij by exact diagonal- 
ization of Eq. ([1} in a finite system with two impurity 
spins either aligned ferromagnetically (FM) or antifer- 
romagnetically (AFM). Then Jij can be expressed as 
the energy difference between the two configurations)^ 
J,^. = [£;(FM) - E{AFM)]/{2S^). We wiU for simplicity 
set S = 1. This solution method is non-perturbative, 
automatically avoids any artificial ultraviolet cut-off de- 
pendencies, and is also capable of generating results for 
any R. However, solving in a finite system creates its 
own problems. We apply periodic boundary conditions 
(PBC) to avoid the effects of edge states but then have 
to deal with the two impurities in one unit cell also in- 
teracting with the impurities in neighboring cells. By 
systemically increasing the "padding" p around the two 
impurities, see Fig. [IJb), we can determine the necessary 
size of the unit cell for converged results for Jij . As seen 
later on in Figs. [3]and[31 p = 2R is in general sufficient. 
We have also ensured convergence with respect to the 
number of fc-points in reciprocal space. 

Apart for studying the RKKY interaction in the bulk 
we have also looked at ZGNRs and strips. Here the 
graphene lattice is terminated along the zigzag direc- 
tion with saturated cr-bonds whereas the p^-orbital on 
the edge atom is unsaturated because of the one miss- 
ing nearest neighbor. We have primarily studied narrow 
symmetric ZGNRs with width W ~ 8/ Via, where a is 
the lattice constant, see Fig. [TJ As is well established, 
the zigzag graphene edge hosts localized states at zero 
energy which significantly changes many of the physical 
properties compared to the bulk.-'^''"^- However, the arm- 
chair edge does not have any such zero energy localized 
states and we find that that for large enough unit cells, 
the results for a ZGNR, where PBC are applied in the 
direction of the ribbon, are the same as those for a strip, 
where instead armchair edges terminate the strip. 

III. RESULTS 

We start with displaying in Fig.[2]the spin polarization, 
— (a|at ~ a.nd s|j, induced into the graphene 

from two impurity spins positioned on-site (a,b) and in 
a plaquette site (c,d), in the FM (a,c) and AFM (b,d) 
configurations, respectively. Below we will discuss each 
of these results. 
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FIG. 2: (Color online) Two impurity spins along the zigzag 
direction positioned on-site on the B-B sublattice (a,b) and 
in the plaquette site (c,d) in the FM (a,c) and AFM (b,d) 
configurations. The impurity spins are marked with white or 
black arrows and the area of the circles on each site is propor- 
tional to the spin polarization, where excess spin-f (spin-^) 
density is red/dark grey (black). 



A. On-site impurities 

An on-site positioning of the impurity spins, where the 
spins sit directly on top of an A or B atom of the graphene 
lattice, has so far been the dominating setup for RKKY 
studies in grapheneJ^^— Since on-site positioning breaks 
the symmetry of the lattice, it is important to distinguish 
between A-A and A-B positioning of the two impurity 
spins. We have studied both of these configurations along 
both the zigzag and armchair directions as well as verified 
our predictions for the asymptotic large- i? behavior for 
several other, chiral, directions. Figures [2][a,b) show typ- 
ical spin polarization patterns of two impurity spins both 
on the B sublattice and along the zigzag direction. We 
clearly see that the spin polarization has different signs 
on the two sublattices close to an impurity and therefore 
one would expect A-A (B-B) impurities to prefer a FM 
coupling whereas AFM coupling should be the case for 
A-B (B-A) impurities. 



Figure |3] shows as a function of the impurity dis- 
tance R for all four different configurations of sublat- 
tice and zigzag /armchair directions for on-site impuri- 
ties. First of all we can directly verify that the RKKY 
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FIG. 3: (Color online) -t/J^Jij for A-A (a,b) and t/JiJij 
for A-B (c,d) on-site positioning of the impurity spins along 
the zigzag (a,c) and armchair (b,d) directions as function of 
impurity distance R. Solid curves represent calculated results 
with padding p = IR — 4i? (black x, blue green *, red o), 
where the lines are only a guide to the eye, whereas the large- 
it dependence in Eqs. (|3I4|) with C — l/(72\/37r) is displayed 
with black Ds. 

coupling is always FM for A-A sites and AFM for A-B 
sites as seen in the sign difference of Jij . This has been 
predicted beforei^ii^ and is a property of the bipartite 
lattice.— Secondly, we conclude that padding p = 2R 
is enough for well converged results. Third, we extract 
the following functional dependence in the large i? = |R| 
limit: 



</ij,A-A(R-) 
Jij,A-B(^) 



Jll + cos(2kD • R) 



C 



t |R|3 

3J| 1 + cos(2ki5 • R + 7r) 



t 



(3) 
(4) 



for R measured in units of the lattice constant a. In 
Fig. [3] these results are plotted as black Ds using the nu- 
merical prefactor C — l/(72\/37r) and, as seen, there is 
essentially a perfect agreement at larger R. In the above 
equations k£) is the reciprocal vector for the Dirac points. 
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i.e. the corners of the Brihouin zone. There are six such 
vectors and for the A-A configuration the resuh is inde- 
pendent on the particular choice of k£) since R is then 
a lattice translation vector, i.e. R = niCi + rt2C2 where 
ni,ri2 are integers. However, for the A-B configuration 
R is not a lattice translation vector and cos(2k£) • R) 
can then depend on the choice of ki). Eq. is only 
correct when k£) • R is maximized, i.e. when k/j is cho- 
sen to be as parallel as possible to R. Eqs. (I3]|4]) are 
very similar to the results derived earlier by Saremiii 
using perturbation theory in a continuum model except 
for the addition of the 7r-phase factor and the need for 
specifying k^ in the A-B result, as just discussed. The 
TT-phase factor is essential for reproducing the numeri- 
cal results for A-B sublattice coupling as these are 180° 
out of phase with the A-A results. These two discrepan- 
cies stem from the fact that in all previous work R has 
improperly been treated as a lattice translation vector 
even for A-B impurities. As seen in our results, choos- 
ing R to be the correct impurity-impurity distance not 
only changes the RKKY coupling at short distances, as 
one might have expected, but also produces the 7r-phase 
shift and a need for an explicit choice of k^i for all R. 
We strongly believe that a proper handling of R in a 
perturbative treatment in the continuum model will also 
include these two additional corrections found in our nu- 
merical results. Our numerical prefactor C — l/(72-\/37r) 
also agrees with the results of Saremiii if one explicitly 
use a factor oi t = 8/3 w 2.67 eV to produce the proper 
energy dimension of their results. We thus conclude that 
our non-perturbative results agree, up to a few small, 
and traceable differences, with the perturbative results 
of Saremi and they firmly establish that the RKKY cou- 
pling is oscillatory for certain R-vectors, although never 
changes sign on the same sublattice, which is contrary to 
some other recent RKKY results.— However, note that 
due to the impurities only appearing at lattice sites, the 
(1 -|-cos(2k£) • R)) -oscillation is in general undersampled. 
For the zigzag direction the period is 3a instead of 3a/4 
whereas for the armchair direction the period is infinitely 
long. Also, for the cases reported in Fig. [3l it is only 
for the zigzag A-B configuration that the RKKY cou- 
pling completely disappears at certain sites. What makes 
the RKKY coupling in graphene unusual is this non-sign 
changing oscillation on the same sublattice as well as the 
decay as compared to Jij oc sin(2/cFi?)/(2fci?i?)^ 
for an ordinary 2D metalj^i^ Finally, we are not only 
able to extract the long-distance behavior but can also 
directly see the deviations from this behavior at short im- 
purity distances. As seen in Fig. [31 the results along the 
zigzag direction are well converged toward the large-i? 
results around R ~ 10a. The exceptions are the results 
for every third lattice site in the A-B configuration which 
are zero in Eq. 21 Along the armchair direction the con- 
vergence is even faster and, except for nearest neighbor 
impurity spins, Eqs. ([3]|4]) give correct results for any R. 
We thus conclude that the large- i? limit is reached for 
surprisingly small impurity distances R. 



B. Plaquette impurities 

For magnetic atoms deposited on graphene, the on- 
site position might not the most energetically favorable 
but instead the atoms can prefer to sit in the middle of 
the hexagon, in the plaquette sitCi ^^'^^ We will here first 
study the simple situation where the impurity spin cou- 
ples incoherently, and with the same coupling constant 
Jk, to all six nearest neighbors in the honeycomb lat- 
tice. Representative spin polarization maps in this case 
are shown in Figs.|2][c,d) which have plaquette impurities 
along the zigzag direction. For this incoherent, symmet- 
ric, coupling to the lattice the large- i? result can in fact 
be directly derived from the on-site results in Eqs. 
by summing the interactions between all combinations 
of nearest neighbors of each impurity spin. Doing so, 
it turns out that in this case the oscillations cancel, the 
AFM coupling prevails, and the asymptotic large- i? re- 
sult is given by 

36 1 

Jij, plaqCR) = C - (5) 

This is the same result as derived by Saremii^ despite 
the additional 7r-phase factor in Eq. (jj]). Figure 21 shows 
our exact diagonalization results together with Eq. ([5l). 
As before, we see that p = 2R is enough to reach con- 
verged numerical results. For plaquette impurities along 
the armchair direction. Fig. Ul^b), there is a systematic 
increase in J^j^piaq for small R compared to the long- 
distance result but the asymptotic behavior is nonethe- 
less reached before R = 10a. For the zigzag direction. 
Fig. HJa), there is some oscillations around the asymp- 
totic value for short R but convergence with respect to 
this value is reached already at i? ~ 4a. 
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FIG. 4: (Color online) — Jij.piaq for impurities in the pla- 
quette site along the zigzag (a) and the armchair (b) direction. 
Same color coding as in Fig. [3l 

While the incoherent case above is straightforward to 
derive from the results of on-site impurities, coherent cou- 
pling, where of a plaquette impurity spin couples to a 
linear coherent combination of the nearest neighbor con- 
duction electrons, is physically more realistic. Saremii^ 
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showed that the same cancelation of the osciUations that 
occur in Eq. (O also makes the 1/B? contribution disap- 
pear altogether for coherent coupling. Since we get the 
same cancellation in the incoherent case when we prop- 
erly treat the impurity distance, we conclude that this 
conclusion is still valid. Coherent plaquette impurities 
thus have a significantly weaker, and thus less relevant, 
RKKY coupling than the other configurations studied 
here. 



C. Impurities in ZGNRs 

Zigzag edges in graphene have proven to be an excit- 
ing new playground for magnetism since this termination 
leads to a localized, flat-band, zero energy edge state 
and is thus extremely prone to spin polarizationii2r— 
This zero energy singularity in the local density of states 
(LDOS) has been claimed to significantly change the 
perturbative RKKY coupling result because of a qual- 
itatively modified behavior of the spin susceptibility at 
the edgeJ^ It is therefore of large interest to further 
study these systems and establish the consequences of 
zigzag edges on the RKKY behavior in an exact, non- 
perturbative setting. We have found that an armchair 
edge does not significantly effect the RKKY coupling and 
thus it is only necessary to study the zigzag edge in order 
to establish the qualitative effect of edges on the RKKY 
coupling. 

Figure [5] shows the spin polarization for two prototyp- 
ical situations in a narrow ZGNR: impurities inside the 
ZGNR (a,b) and along the edge (c,d). The RKKY cou- 
pling — Jij for the configurations in Fig. [S] is shown in 
Fig. [D^a) for the parameters t — Jk ^ I- While this is 
a rather unphysically high value for Jk it helps display- 
ing the essential features in a numerically accessible R- 
range. We have for comparison also included the results 
for Jk = t/100, and, as discussed below, the same physi- 
cal behavior is governing both J^-values, although all sig- 
nificant features get extended over a longer i?-range for 
smaller Jk, making a complete numerical study harder. 
For reference in Fig. [51 the asymptotic zigzag result in the 
bulk is shown in black and we directly see that impuri- 
ties in narrow ZGNRs behave significantly different. The 
coupling between edge impurities (red, o) is larger than 
in the bulk for small distances and also displays some os- 
cillations in that regime, but then decays exponentially 
for larger R which is in sharp contrast to the power-law 
decay in the bulk. Displayed is also Jij for the RKKY 
coupling between opposite edges, the opposite sign be- 
ing a consequence of the two edges belonging to different 
sublattices. As seen, after an initial short distance, where 
the width of the ZGNR is important, same edge and op- 
posite edge impurities couple with equal strength. This 
is another difference compared to the bulk results where 
the AFM coupling is 3 times larger than the FM cou- 
pling. We thus draw the conclusion that the edges are 
dominating the response for edge impurities. This could 




FIG. 5: (Color online) Two on-site impurity spins inside a 
narrow ZGNR (a,b) and on the edge of the same ZGNR (c,d) 
in the FM (a,c) and AFM (b,d) configurations, respectively. 
Same color coding as in Fig. [21 



have been anticipated already from Figs.[5ljc,d) where the 
spin polarization is seen to be large only in the vicinity of 
the impurity and only along the edge, it does not spread 
significantly into the bulk. The essential features of the 
RKKY coupling along the edge can be understood from 
the following argument: A magnetic impurity will always 
force a spin polarization of the graphene in its vicinity. 
For a zigzag edge impurity this spin polarization can triv- 
ially, and essentially without energy penalty, be achieved 
by polarizing the localized edge state at zero energy with- 
out much polarization of the surrounding bulk. Conse- 
quently, away from the impurity, the edge state will also 
quickly become unpolarized, much more quickly than the 
bulk can lose its spin polarization, and thus the RKKY 
coupling between two edge impurity spins decays much 
faster than in the bulk. With this argument one would 
expect the total amount of polarization in the B sub- 
lattice for the FM configuration, ^ |sb|, to be large for 
small R, as then the two impurities interact very strongly 
with each other, but then rapidly decay with R until its 
flattens out to a value equal to the sum of the total spin 
polarization of two uncoupled edge impurity spins. This 
behavior is confirmed in Fig. where the asymptotic 
value is shown to be reached already around R — 7a 
for Jk — t. This should be compared to the evolution 
of the spin polarization for the equivalent configuration 
in the bulk where the total spin polarization is almost 
constant, as displayed by the black curve in Fig. [HJb). 
The exponential decay of Jij is in sharp contradiction to 
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FIG. 6: (Color online) (a) —Jij for impurities along the zigzag 
direction on the zigzag edge of a narrow ZGNR (red, o), inside 
the same ZGNR (green, and the asymptotic behavior in 
the bulk (black, x) as function of impurity distance R. In 
(blue, +) is Jij for impurities on opposite edges of the ZGNR. 
Here t = Jk = 1 expect the gray line where t = 1, Jt = t/100. 
(b) The total spin polarization in sublattice B, for the 

same situations as in (a) in the FM impurity configuration. 
A large enough unit cell was used to ensure convergence of 
^ |ss| in the 7?-range displayed. 



earlier calculations on the same vifidtli ZGNR by Bunder 
et al}^ who used an analytical approach to the tight- 
binding structure of graphene to calculate the spin sus- 
ceptibility and thus obtained perturbative results for the 
RKKY coupling. Both methods predict an equivalence 
between same and opposite edge impurity spins but Bun- 
der et al. report an almost linear decay with distance 
for small distances followed by sign oscillations in the 
RKKY coupling. Despite a thorough investigation of the 
RKKY coupling for R < 130a, we were not able to detect 
any deviations from the exponential decay, and thus no 
sign changes, within the numerical accuracy which was 
a factor of 10"" of the RKKY couphng at i? ~ 0. We 
thus conclude that when studying these edge states, a 
non-perturbative method appears to be essential when 
calculating the RKKY coupling. 

In the bulk Jij cx and X^ksl oa Jk but for edge 
impurities these simple relations do not longer hold. 
For all Jk > t/10 we have found that the asymptotic 
large- i? induced spin polarization from edge impurities is 
^ |sb| « 1.1 and even for Jk = t the asymptotic value 
of ^ I SB I is only slightly higher due to a finite amount 
of induced polarization in the nearby bulk, see Fig.lSJb). 
Thus even in the limit of vanishing Jk, edge impurities 
are going to elicit a finite spin polarization response of 
the ZGNR. This is again due to the extreme easiness of 
polarizing the zero energy edge state. However, when 
Jk decreases the spin polarization per edge site naturally 
goes down so the polarization now instead have to be 
spread over more edge sites. This has an interesting con- 
sequence for the RKKY coupling: When Jk decreases, 
the more elongated polarization response causes the os- 
cillations in the RKKY coupling for small R to be spread 



out over a longer i?-range. Eventually, however, an expo- 
nential decay is achieved for large enough R for any Jk, 
as also seen in the gray curve in FigHl^a) for Jk — t/100. 
But note that the exponential decay rate also becomes 
smaller when the spin polarization gets more elongated 
along the edge, thus making the RKKY coupling at large 
R larger for decreasing Jk- As a direct consequence, the 
RKKY coupling between edge impurities is going to be 
larger than the coupling between bulk impurities over a 
larger i?-distance the smaller the Jk- Of course in the ex- 
treme large- i? limit the exponential decay is always going 
to make the edge impurity RKKY coupling smaller than 
the bulk coupling. It is worth pointing out before leaving 
the treatment of same edge impurities that for all of the 
effects described above, the width of the ribbon is not 
essential and we thus expect the same behavior even for 
much wider ribbons. 

Physically in between zigzag edge and bulk impuri- 
ties we find impurities positioned inside a narrow ZGNR. 
The spin polarization for this situation is displayed in 
Fig. [Sl[a,b). Here the bulk is polarized in the direction 
parallel to the ZGNR but this bulk polarization also in- 
duces an edge polarization. Thus the absolute value of 
the polarization is here very large and it grows with the 
distance R since all the bulk between the two impurities 
is at least lightly polarized and any amount of bulk po- 
larization is going to elicit a large polarization of the edge 
state. This behavior is clearly seen in Fig.lGjb) where the 
total spin polarization increases linearly with R. How- 
ever, the edge states are also very easily unpolarized, 
as discussed above, so the effect on the RKKY coupling 
from the edges is not expected to be very large compared 
to the bulk contribution. This prediction is confirmed in 
Fig. [nja) where we see that, while the coupling is some- 
what larger than in the bulk and lacks oscillations, it 
decays with approximately the same power-law as in the 
bulk. Also, Jij oc and X^l'^sl Jk for these impu- 
rities, which further establish the bulk-like behavior of 
impurities inside a narrow ZGNR. We anticipate that for 
wider ZGNRs all bulk properties are fully recovered. 



IV. DISCUSSION 

Our results above are obtained without any approxi- 
mations except the finite size of the system, assuming a 
non-interacting picture for the electrons in graphene, and 
using a nearest-neighbor hopping band structure. The 
first approximation we have taken care to handle system- 
atically and as long as the padding p around the impu- 
rities are twice or larger than the impurity-impurity dis- 
tance R in each direction, the results are well converged. 
Also, since the asymptotic behavior is reached for rela- 
tively small R we have been able to extract the large R 
limit to make a comparison with earlier, partly disagree- 
ing, resultai^^— obtained using a perturbative approach 
within the continuum field-theory model. Eqs. (|31)-(|31) are 
our results for large R for bulk impurities and these are 
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closely related to one of the results in the literature's al- 
though a TT-phase factor and the need for explicitly choos- 
ing fc^) for the A-B configuration are new findings. This 
difference is most likely due to a previous improper treat- 
ment of the impurity distance for impurities on differ- 
ent sublattices. Thus our results establish both that the 
usual perturbative treatment of RKKY coupling is appro- 
priate in the bulk and, that, while the RKKY coupling 
in graphene does not undergo sign changes with distance 
on the same sublattice, it still has a (1 -I- cos(2k£) • R))- 
oscillation, a fact that has only been pointed out in one of 
the previous works.^^ Despite these oscillations the facts 
that the AFM A-B coupling is 3 times larger than the 
FM A-A coupling and that the coupling is stronger at 
short distances result in a strong tendency towards AFM 
order for any random configuration of impurities. Thus 
these new results still support the conclusion of mag- 
netic defects in graphene creating a dilute antiferromag- 
net at low enough temperatures.^^ The moments will here 
be oriented in opposite directions on the two sublattices 
with the total magnetic moment equal to zero. Since the 
RKKY coupling is always AFM (incoherent coupling) or 
very small (coherent coupling) for plaquette impurities 
a dilute AFM state could also be present for plaquette 
impurities. We have also studied ZGNRs where we have 
shown that the zero energy edge state present on the 
zigzag graphene edge significantly modifies the coupling 
between impurity spins. We show that defects along the 
edges couple very strongly at short distances but that 
the coupling finally decays exponentially with distance 
in contrast to the power-law decay in the bulk. The 
exponential decay is a consequence of the extreme eas- 
iness of polarizing and unpolarizing the localized zero 
energy state. This result disagrees with earlier perturba- 
tive results'^ where they found the decay to be almost 
linear for small R followed by oscillatory sign changes. 
So while the standard perturbative RKKY treatment is 
approximately valid in the bulk we have here shown that 
for ZGNR edge impurities a numerical non-perturbative 
treatment is essential. 

The approximation of ignoring electron-electron inter- 
action in graphene is on the other hand harder to moti- 
vate. The non-interacting picture has been predominant 
in the study of RKKY couplingi^"— as well as in other 
theoretical work on magnetic adatoms in graphenci^"— 
However, there seem to be growing evidence for the 
importance of electron-electron interactions in graphene 
with theoretical results pointing to graphene being close 
to a Mott insulator state<^""— Electron interactions could 
thus potentially significantly modify the magnetic prop- 
erties of graphene and therefore also the RKKY coupling. 
One such example would be if the insulating state is a 
Neel phase*^ While a comprehensive treatment of elec- 
tron interactions in graphene is extremely hard, a density 
functional theory (DFT) calculation of a carefully chosen 
system should be able to offer insight into the impor- 
tance of interactions for the RKKY coupling. However, 
not only will the magnetic moment and its coupling to 



graphene be more complicated in a real system than in 
the idealized Eq. ([T]) but also a large unit cell is needed, 
making the DFT calculation highly computationally in- 
tensive. Nonetheless, some DFT results exist for short 
distances^° and line defects, both pointing to a slower 
power- law decay than R~^. A detailed comparison of 
such results with our short distance RKKY coupling may 
offer insights in the applicability of the non-interacting 
approximation and is the goal of a future study. In fact, 
initial data points to the possibility of electron-electron 
interactions producing strong enough corrections to the 
non-interacting results that it is reasonable to then ig- 
nore any corrections to the nearest neighbor hopping 
band structure in comparison. In a ZGNR the effect of 
electron-electron interaction is possibly larger than in the 
bulk as here any electron-electron interaction will drive 
a spontaneous spin polarization of the edges.—"— Thus, 
in a real ZGNR the edge is already polarized and an im- 
purity spin can then not as easily polarize the edge as 
found in the non-interacting picture, especially since Jk 
is usually a small parameter. A recent DFT study^ has 
shown on an almost exponentially decaying RKKY cou- 
pling for nearby impurities inside a narrow ZGNR which 
points to the importance of electron-electron interactions 
even inside ZGNRs. This would be in contrast to the 
non-interacting picture where impurities inside a ZGNR 
behave essentially bulk-like. Detailed results on the influ- 
ence of electron-electron interactions for ribbons is also 
a subject of the future study. 

We have here also explicitly ignored the issue of how an 
external magnetic moment is created in graphene, but for 
any experimental verification or use of our results such 
considerations have to be taken into account. The sim- 
plest implementation is probably to deposit a magnetic 
adatom on top of graphene, such as Co or Fe^^^ We 
also expect the long distance behavior to be qualitatively 
similar for substitutional magnetic atoms, of which at 
least Co has been shown to be magnetic in graphene.— 
Single-atom vacancies and hydrogen chemisorption de- 
fects have also both been shown to possess a magnetic 
moment in graphene^"~ and should behave similarly to 
substitutional magnetic atoms. In fact, for vacancies, 
Lieb's theorem*^ gives the same AFM/FM state as found 
here. It has even been proposed that the FM state found 
in proton-radiated graphite^^^ is due to excess vacancy 
creation on one sublattice, which is in agreement with our 
findingSi-t^ At very short impurity distances, however, an- 
nihilation of vacancies as well as direct exchange between 
magnetic moments can also be of importance. In aggre- 
gate, there exist multiple possibilities for experimentally 
studying the RKKY coupling in graphene. Moreover, the 
critical coupling for the Kondo effect is likely too high in 
undoped graphen o^^i^^'^°i^° and therefore the Kondo ef- 
fect should not compete with the RKKY coupling in any 
of these systems. 

In summary we have studied the RKKY coupling be- 
tween two impurity spins on graphene for any impu- 
rity distance R using exact diagonalization. Our results 
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largely agree with an earlier perturbative result in the 
large R limit where oc ±(H-cos(2k£) •R))/|Rp, with 
A-A sublattice arrangement FM and A-B sublattice ar- 
rangement AFM, although an additional 7r-phase shift for 
the A-B sublattice arrangement is a new finding. The 
large R limit is reached within a few lattice unit con- 
stants and the deviations are in general not large even 
for small R in the bulk. We have also studied the ef- 
fect of zigzag edges and found that impurities along this 
edge display an exponentially decaying coupling at large 
R due to the easiness of polarization of the zero energy 
edge, state in an non-interacting picture. This result for 
edge impurities is however in contrast with earlier pertur- 



bative results/^ pointing to the importance of an exact, 
non-perturbative, treatment of the RKKY coupling in 
ZGNRs. Impurities away from the edge, even in narrow 
ZGNRs, regain most of the bulk characteristics. 
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